Solitons and their ghosts in T'T-symmetric systems with defocusing nonlinearities 
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We examine a prototypical nonlinear Schrodinger model bearing a defocusing nonlinearity and Parity-Time 
(VT) symmetry. For such a model, the solutions can be identified numerically and characterized in the pertur- 
bative limit of small gain/loss. There we find two fundamental phenomena. First, the dark solitons that persist 
in the presence of the PT-symmetric potential are destabilized via a symmetry breaking (pitchfork) bifurca- 
tion. Second, the ground state and the dark soliton die hand-in-hand in a saddle-center bifurcation (a nonlinear 
analogue of the PT-phase transition) at a second critical value of the gain/loss parameter. The daughter states 
arising from the pitchfork are identified as "ghost states", which are not exact solutions of the original system, 
yet they play a critical role in the system's dynamics. A similar phenomenology is also pairwise identified 
for higher excited states, with e.g. the two-soliton structure bearing similar characteristics to the zero-soliton 
one, and the three- soliton state having the same pitchfork destabilization mechanism and saddle-center collision 
(in this case with the two-soliton) as the one-dark soliton. All of the above notions are generalized in two- 
dimensional settings for vortices, where the topological charge enforces the destabilization of a two-vortex state 
and the collision of a no-vortex state with a two-vortex one, of a one- vortex state with a three- vortex one, and 
so on. The dynamical manifestation of the instabilities mentioned above is examined through direct numerical 
simulations. 



I. INTRODUCTION. 

Over the past decade, and since they were originally pro- 
posed by Bender and co-workers (H Q], systems character- 
ized by PT-symmetric Hamiltonians have become a subject 
of intense research efforts. The interest on these systems 
arises from their fundamental property to exhibit real spec- 
tra, while non-Hermitian, thus providing an intriguing alter- 
native to standard Hermitian quantum mechanics. In the case 
of a standard Schrodinger-type Hamiltonian, with a generally 
complex potential U, the VT symmetry dictates that the po- 
tential satisfies the condition U(x) = U*(—x), where (•)* 
stands for complex conjugation. 

An important development in the study of such models was 
the work of Christodoulides and co-workers SH], who pro- 
posed and studied (both theoretically and experimentally) op- 
tical systems featuring the VT symmetry. In the optics con- 
text, however, another key element that comes into play in the 
physics of such systems is nonlinearity. Therefore, the con- 
siderations of Refs. ||3l|4|] extended from bright and gap soli- 
tons to linear (Floquet-Bloch) eigenmodes in periodic poten- 
tials, examining how these coherent structures are affected by 
the genuinely complex, yet PT-symmetric potentials. More 
recently, experimental results were reported both in nonlin- 
ear optical systems (HEI and electronic analogs thereof 0]. 
These, in turn, have triggered a wide range of theoretical stud- 
ies on nonlinear lattices with either linear 1181-12 ill or nonlin- 
ear ll22l - [24ll PT-symmetric potentials and, more recently, on 
harmonic PT-symmetric potentials l25ll . 

In the above works, numerous features extending from 
bright solitons to defect modes, and from gap solitons to PT- 
lattices have been examined. Nevertheless, the consideration 
of defocusing nonlinearities, and especially of dark solitons 
has received limited attention; see, e.g., Refs. |26|,|27|]. An- 
other theme that despite its considerable relevance has also 



been considered only by a few works is that of higher dimen- 
sionality. For the latter, the focusing nonlinearity case JH 0] 
has been examined, especially so in the context of lattice po- 
tentials. Our aim herein is to provide a systematic analysis 
of PT-symmetric Hamiltonians exhibiting defocusing non- 
linearities. Particular, and building on the earlier work of 
Ref. [27], we give a detailed account of the existence, sta- 
bility and dynamical properties of the ground state and first 
few excited states namely dark solitary waves and vortex two- 
dimensional generalizations thereof. This is done specifically 
in a prototypical context of the nonlinear Schrodinger type 
which can, in principle, be relevant both in the case of Bose- 
Einstein condensates l28h . as well as in that of nonlinear op- 
tics J29D; see the detailed justification of the real and imagi- 
nary parts of the potential below. 

Our main findings and their presentation is structured as 
follows. 

• In Sec. II, we put forth the general model and consider 
its physical relevance. 

• In Sees. III-IV, we turn to generic potentials and the 
features of their ground and excited states in one- 
dimension. We find that odd excited states (1 -soliton, 
3-soliton, 5-soliton, etc.) become subject to a sym- 
metry breaking bifurcation. The pitchfork nature of 
this event is rationalized through the introduction of so- 
called ghost states which arise from it and which are 
exact solutions of the steady state but remarkably not 
ones such of the original full problem. Nevertheless, 
direct numerical simulations clearly identify the ensu- 
ing symmetry breaking and manifest the dynamical role 
of the ghost waveforms. One more important feature 
identified is the nonlinear analogue of the PT-phase 
transition. In particular, between the (unstable, for suf- 
ficiently strong gain/loss) saddles of the odd excited 
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states and the centers of the even states (ground state, 2- 
soliton, 4-soliton, etc.), there is a pairwise collision and 
disappearance (blue- sky or saddle-center) event. This 
is strongly reminiscent of the corresponding transition 
of linear eigenstates of the Hamiltonian originally re- 
ported in Ref [1]. Direct numerical simulations are em- 
ployed in order to identify the evolution dynamics of 
the unstable solitary waves states and also their dynam- 
ics past the PT-phase transition point in Sec. V. The 
ghost states are then examined separately in their own 
right in Sec. VI. 

• In Sec. VII, we illustrate how each of the above parts of 
the picture is generalized in the two-dimensional vari- 
ant of the relevant model. There, the charge of the vor- 
tex states imposes topological constraints enforcing that 
the ground state may only collide and disappear (in a 
nonlinear PT-phase transition) with the two-vortex (of 
opposite charge, namely dipole) state. Similarly, the 
single vortex and triple (of alternating charge) suffer 
a saddle-center bifurcation and so on. Prior to these 
events, a destabilization of the two- vortex (and three- 
vortex etc.) states arises through a pitchfork event creat- 
ing ghost vortex states. The latter are illustrated dynam- 
ically as well and the structural analogies of the one- 
and two-dimensional settings are explored both in the 
statics and in the dynamics. Lastly, Sec. VIII presents 
a brief summary of our conclusions and a number of 
potential directions for future study. 



II. THE MODEL AND SOME ANALYTICAL INSIGHTS 

Our model, which can be equally applied to a variety of 
one-dimensional and even higher dimensional systems is a 
nonlinear Schrodinger (NLS) equation incorporating a com- 
plex potential. This equation, which finds applications in the 
contexts of nonlinear optics [29] and in the physics of atomic 
Bose-Einstein condensates 12811 . is expressed in the following 
dimensionless form: 

id t u = ~\ d l u + \ u ? u + [V(x) + iW(x)]u, (1) 

where u is a complex field, denoting the electric field enve- 
lope in the context of optics (or the macroscopic wavefunc- 
tion in BECs), t denotes the propagation distance (or time in 
BECs), x is the transverse direction, while V(x) and W(x) 
denote, respectively, the real and imaginary parts of the ex- 
ternal potential. For a VT- symmetric Hamiltonian, V(x) and 
W{x) must be, respectively, an even and an odd function of 
x, namely: 

V(x) = V(-x), 

(2) 

W{x) = -W(-x). 

Physically speaking, in the context of optics, V(x) and W(x) 
represent, respectively, the spatial profiles of the real and 



imaginary parts of the refractive index. The above require- 
ments for the parities of V(x) and W(x) can be met in a case 
where V(x) has, e.g., a parabolic profile and W(x) has an 
anti- symmetric profile (amounting to equal and opposite gain 
and loss), as in the experiment of Ref. [6]. In the context of 
BECs, V(x) represents the external trap (necessary to confine 
the atoms 1 28]) and W(x) accounts for a mechanism for in- 
jecting and removing particles in equal rates. The requirement 
for the parity of the trap V(x) can easily be met in the case 
of e.g., the usual parabolic potential (representing a magnetic 
trap) or a double-well potential (representing a combination 
of a magnetic trap with a suitable optical lattice) [28]. On 
the other hand, W{x) is odd if an equal number of atoms is 
injected and removed from spatial regions symmetrically lo- 
cated around the trap center. Such a setting was originally 
proposed in Ref. [30] for a BEC confined in a double-well po- 
tential. Furthermore, passive parity-time symmetric analogs 
of the double well system (featuring only loss in one well and 
no gain) have also been proposed in the context of the so- 
called open Bose-Hubbard dimer [ 3 1 ] ; see also Refs. iHd. 
More recently, a very large volume of activity has focused on 
such double- well potentials in a balanced gain-loss (i.e., VT- 
symmetric) form; see, for instance, Refs. HH - HH] . In such 
contexts, some of the notions presented below, such as the 
ghost states and their emergence from suitable bifurcations 
and dynamical relevance, can be both numerically manifested, 
as well as analytically demonstrated. 

We seek standing wave solutions of Eq. (Q]) in the form 
u = p(x) ex.p[i(j)(x) —ifit], where the real functions p(x) and 
(j>{x), and the real constant p represent, respectively, the am- 
plitude, phase and propagation constant (in optics) or chemi- 
cal potential (in BECs). Substituting this ansatz into Eq. (Q]), 
and separating real and imaginary parts, we obtain the follow- 
ing coupled boundary-value problems (BVPs): 

MP = -7>(pxx ~ p4>l) + P 3 + V(x)p (3) 

2W(x)p 2 = (p 2 $ x ) x , (4) 

where subscripts denote partial derivatives. For generic poten- 
tials V(x) and W(x), the relevant PT-symmetric problem in- 
volves solving the BVPs for p and <fi. Notice the critical role of 
the imaginary part of the potential in dictating the phase pro- 
file (j)(x). We will chiefly focus on the case of a real parabolic 
potential, 

V(x) = \tfx\ (5) 

with strength Q, modeling the transverse distribution of the 
refractive index (or the external trap in BECs) as mentioned 
above, while the imaginary part W(x) will be considered to be 
an odd, localized function of space, of spatial width <0 _1 . 
Such a form of W(x) is consistent with the experimental work 
of Ref. |6|] in a nonlinear optics setup and could be relevant 
to an effective description of atom loss/gain mechanisms in 
trapped BECs. Our analysis will be general (independent of 
the particular form of W(x)) 9 and we will showcase our re- 
sults in the case of the following prototypical example: 

W(x) = sxexp(-x 2 /2), (6) 
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FIG. 1: (Color online) The density [solid (blue) line] and phase [solid 
(red) line] of the numerically obtained TF background compared to 
the prediction of Eqs. (UTT) [dashed (green) lines]. The parameters 
values are: \i — 3, f2 = 0.1 and e = 0.3. 



where e is a parameter setting the magnitude of the imaginary 
potential (a generalization of this model in two-dimensions 
will be studied in Sec. I VIII) . We should mention that our an- 
alytical approximations (see below) can also be applied for 
other choices, e.g., when W(x) takes the form W(x) = 
esech 2 (x) tanh(x); in such cases, we have checked that our 
results remain qualitatively similar to the ones that we will 
present below. 

We conclude this section by noting the following. Below 
we will study the ground state and excited states of the system 
(in the form of dark solitons). In that regard, it is relevant 
to consider the evolution of the physically relevant quantity 
N = J \u\ 2 dx, which represents power in optics or number 
of atoms in BECs. Employing Eq. (Q]), it is straightforward to 
find that dN/ dt is governed by the equation: 

— = 2 J \u\ 2 W(x)dx. (7) 

Thus, since W{x) is odd, it is obvious that the power is con- 
served as long as the square modulus profile of the ground 
state or of the excited states remains even. Below, we will 
show that for genuinely stationary states this is the case, in- 
deed and we will examine the important consequences of 
Eq. © on the bifurcations and dynamics of the system at 
hand. 



III. GROUND STATE AND SINGLE DARK SOLITON 

A. Ground state 

First, we will study the existence and stability of the most 
fundamental state of the system, namely the ground state. The 
latter is sought as a stationary solution of Eq. (Q]) in the form 
u = Ub{x) exp(— ifit) (where fi is the propagation constant or 
the chemical potential in BECs), with the background field u b 



obeying the equation: 

- \ a l u *> + \u b \ 2 u b + [V(x) + iW(x)]u h - fiu h = 0. (8) 

For a sufficiently small imaginary potential, W(x) = eW(x) 
[with max{|W"(x)|} = 0(1)], where e <C 1, and when 
the inverse width of V(x) is sufficiently large so that 
Q ~ £, we may find an approximate solution of Eq. (Q]) in 
the Thomas-Fermi (TF) limit (HHH)]. This solution is of the 
form: 

u h (x) = [y/Ji + / (x)} exp[i</>(x)], (9) 

where the amplitude and phase f(x) and <j)(x) (considered to 
be small, of order e 2 and e, respectively) are given by: 

/(*) = -^=( F + 2W2 )' ( 10 ) 

(j){x) = 2 Jwdx. (11) 

where W = J Wdx, and we have neglected terms of order 
0(e 3 ) [notice that the integral in Eq. (TTTb is an indefinite one]. 
Contrary to the conservative case (e = 0) [28], this TF back- 
ground is characterized by a density dip located at the cen- 
ter, x = 0, and a nontrivial (tanh- shaped) phase distribution. 
Both features are shown in Fig. [TJ where the analytical re- 
sult (dashed lines) is compared with the numerical one (solid 
lines); the agreement between the two is excellent. Note that 
we have checked that the above solution stays close to the 
numerically found ground state of the system up to the order 

o(e 2 ) a 

The evolution of the density of the TF background indi- 
cates that the ground state is stable (see more details in the 
analysis below). Importantly, a linear stability — Bogoliubov- 
de Gennes (BdG) — analysis (see, e.g., Ref. IK38I1 ) justifies the 
above argument, showing that the background u b (x) is indeed 
stable against small perturbations. 

B. Single dark soli ton 

Apart from the ground state, excited states of the system 
— in the form of stationary dark solitons — can also be found 
numerically, by means of a fixed point algorithm (Newton's 
method). A pertinent example of the form of a single dark 
soliton is shown in the top panel of Fig. [2 

To analyze the dynamics of such a single dark soliton, de- 
scribed by function v(x, t), on top of the TF background, we 
introduce the product ansatz: u = Ub(x)v(x,t) into Eq. (Q]) 
(i.e., we decompose the solution into a soliton and its back- 
ground). Then, we employ the expressions (f9t— (TTTb for u b , 
and simplify the resulting equation for v(x, t) by keeping only 
leading-order terms, up to order G(e 2 ) [recall that the function 
f(x) is of order G(e 2 )]. To this end, using the scale trans- 
formations t fit and x —> y/Jix, we obtain the following 
perturbed NLS equation: 

idtv + -d 2 x v + - M 2 ) = P(v), (12) 
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FIG. 2: (Color online) Top panel: The density [solid (blue) line] 
and phase [dashed (red) line] of the first excited (single dark soliton) 
state for e = 0.3. Middle and bottom panels: The linear spectrum of 
the single dark soliton branch: middle (bottom) panel shows the real 
(imaginary) part of the lowest eigenfrequencies cj, with respect to the 
amplitude e. The dashed (red) line depicts the mode u) a , which co- 
incides with the anomalous mode for e — 0; while the dotted (black) 
line is the analytical result of Eq. (EDI . The vertical line shows the 
point eil\ where the bifurcation emerges. The parameters used are: 
[i = 3 and Q = 0.1. 



lution equations for xo(t) and cp(t) read: 



dx . , . 
— = sm,(t) 

dt 2 



(15) 



sech 4 (0 [tanh(0>V 2 + WW] dx, 

(16) 



where we have assumed almost black solitons with suffi- 
ciently small phase angles. Thus, generally, for a given imag- 
inary potential W(x), and by calculating the above integral, 
we can derive an equation of motion for the dark soliton cen- 
ter xq in the form: 



d 2 xo 



dv e 



eff 



dx 



(17) 



where V G s, is the effective potential felt by the single soliton. 

In the case where the imaginary part of the potential takes 
the form of Eq. ©, we numerically calculate the integral in 
Eq. (IT6l> and find that the effective potential can be approxi- 
mated as: 



V e fs(xo) 



71 



£ 2 Ci 

4c 2 



sech 4 (c2£o). (l^) 



where c\ w 1.9 and C2 ~ 0.6. Notice that the effective poten- 
tial contains essentially two contributions, the first from the 
external parabolic trap (as is the case in the context of BECs 
l22l[38ll40]a), and the second from the presence of the imag- 
inary potential W(x). 

To examine the stability of the equilibrium at xq = (i.e., 
the stability of a quiescent black soliton located at the trap 
center), we Taylor expand the effective potential of Eq. (TT8b 
around xq = 0, and find — to leading order — the following 
equation of motion for the soliton center 112711 : 



where the functional perturbation P(v) is given by 



P(V) = fJL- 



(1-M 2 )*; (F + 2W 2 ) 



+v x ( -V x - 2(W - i)W 



Applying the perturbation theory for dark solitons 13 8L 14111. 
we seek a solution of Eq. (TT"2l) in the form of the dark soliton 
of the unperturbed system (P(v) = 0): 



v(x, t) = cos (f(t) tanh £ + i sin (p(t). 



(14) 



Here, £ = cos (p(t) [x — xo(t)], while the unknown, slowly- 
varying functions ip{t) and xo(t) represent, respectively, the 
phase (\(p\ < it/ 2) and center of the soliton. In the framework 
of the adiabatic approximation, the perturbation-induced evo- 



d 2 x 



2 

-0J QSC Xq, 



-p ) - l.Ue z 



(19) 
(20) 



Equation (|20l) implies that if the amplitude e of W(x) is less 

than a critical value « 0.66 Q, the soliton performs oscil- 
(13) lations in the complex potential with frequency u osc ; on the 
other hand, if e > eil^ the soliton becomes unstable. In fact, 
the prediction has the quintessential characteristics of a pitch- 
fork bifurcation. For e < , the relevant potential is monos- 
table at xo =0, while beyond the critical point, two symmet- 
ric minima arise (in this effective potential picture) and the 
center position becomes the saddle point that separates them. 

The above prediction about the change of stability of the 
critical point at xq = has been confirmed numerically, both 
by means of direct simulations and by employing a BdG anal- 
ysis. The latter reveals that the considered stationary dark soli- 
ton is characterized by a mode u a (initially coinciding with 
the anomalous mode (38l [40] for e = 0), which is real for 



e < e, 



(i) 



(in this case, u a = cj OS c), and it becomes imaginary 
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for e > Scl\ thus signaling the onset of the spontaneous sym- 
metry breaking (SSB) instability of the dark soliton (which 
should be expected to displace the soliton from the trap cen- 
ter). 



The dependence of u a on the amplitude e of the imaginary 
potential W, as found by the BdG analysis, is illustrated in 
the middle and bottom panels of Fig. [2] The mode, corre- 
sponding to the oscillatory motion of the dark soliton, is the 
first non-zero mode in the real part of the spectrum, indicated 
with a dashed (red) line — cf. middle panel of Fig. [21 We ob- 
serve that it initially moves towards the spectral plane origin, 
and past the critical point, e& (cf. vertical line), uj a becomes 
imaginary, and thus the soliton becomes unstable. The corre- 
sponding pair of imaginary eigenfrequencies is shown in the 
bottom panel of Fig. [2] The solid (black) line in both panels, 
shows the analytical result of Eq. (l20b . An excellent agree- 
ment between the analytical prediction and the BdG numeri- 
cal result is observed even beyond the SSB bifurcation point, 
while for larger values of the parameter e, perturbation theory 
fails and, as expected, the agreement becomes worse. 



Importantly, for large values of the parameter e, the imagi- 
nary eigenfrequencies start moving towards the spectral plane 

(2) 

origin and collide with it at a second critical point, Scr , af- 
ter which the branch is terminated. To better understand how 
the branch ceases to exist, we first observe (bottom panel of 
Fig-[B that the density profile of the soliton becomes increas- 
ingly shallower (i.e., more "grey") as e grows and the sec- 
ond critical point is approached. This is due to the devel- 
opment of an increasingly strong even imaginary part of the 
solution. Furthermore, the stable background solution Ub(x) 
[cf. Eqs. (ITTb and top panel of Fig. [U develops an odd imagi- 
nary part resembling a (progressively darker) grey soliton. Fi- 

(2) 

nally, at e = Scr , the profiles of these modes become identi- 
cal and disappear in a blue- sky bifurcation through their col- 
lision. In this saddle-center bifurcation, the dark soliton plays 
the role of the saddle, while the ground state of the system is 
the center. This observation leads to the conclusion that there 
is a phase transition where the system loses its PT-symmetry, 
due to the nonlinear problem eigenvalues ceasing to be real, in 
direct analogy with the linear VT phase transition of Ref. Q)]. 
This will be discussed in more detail below. 



However, we should highlight here that one more nontriv- 
ial question persists in the context of the single dark soliton. 
Past the critical point , it is expected that the instability of 
the soliton at the center will provide us with "daughter states" 
which are spontaneously symmetry broken dark solitons cen- 
tered at a positive or at a negative value of xo . For real chem- 
ical potentials/ propagation constants, such states can, never- 
theless, not be identified. This appears to be inconsistent with 
the pitchfork nature of the bifurcation and constitutes a point 
to which we will return after we visit the dark soliton evolution 
to appreciate the dynamical implications of this instability. 




FIG. 3: (Color online) Top panels: the density (left) and phase 
(right) profiles of a stationary two-dark- soliton state. Middle panel: 
the dependence of the equilibrium configuration of the two- soliton 
solution as a function of e\ dashed (white) lines show the theoreti- 
cal predictions for the position of the centers of the waves, cf. fixed 
points of Eqs. j23t-j24l. Bottom panels: the linear spectrum of the 
two soliton branch, namely the real (top) and imaginary (bottom) 
parts of the eigenfrequencies as functions of e. Solid red and green 
lines, correspond to the analytical result of Eq. (l29l . The value of 
these modes at e = 0, corresponds to the first and second anomalous 
modes of the Hamiltonian case. 



IV. MULTIPLE DARK SOLITON STATES AND 
NONLINEAR VT PHASE TRANSITIONS 

In this section, we will study higher excited states, in the 
form of multiple dark solitons. Such states are known to exist 
in the presence of V(x), and have been studied extensively in 
the context of BECs [38, 40151. Here, we will consider them 
in the presence of the imaginary part of the potential W(x). 
Particularly, we will focus on the two- and three-dark- soliton 
states; qualitatively similar results can also be obtained for 
higher excited states. 



A. Two-dark-soliton state 

The density and phase of a two-dark- soliton- state are shown 
in the top panel of Fig. [3] Similarly to the ground state branch 
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(see Fig.Q]), we observe that a density dip is formed at the cen- 
ter, which becomes deeper as e is increased. Additionally, we 
observe that the two solitons are located away from the cen- 
ter, at some finite equilibrium distance ±X. This equilibrium 
corresponds to the situation where the repulsive force between 
the solitons and the "expulsive" (with respect to the trap cen- 
ter) effect, for sufficiently large e, of the imaginary part of 
the potential W(x) analyzed in the last section is counter bal- 
anced by the confining nature of the real part of the potential 

v{x) SH. 

Below we will study analytically how the presence of the 
imaginary potential W(x) modifies the equilibrium distance 
X characterizing the stationary two-soliton-state, and also 
study small-amplitude oscillations of the solitons around their 
fixed points. To obtain an approximate equation of motion 
[analogous to Eq. (l20l) 1 for the two solitons centers, we as- 
sume that each of the two solitons feels the effective external 
potential of Eq. (IT8T) and, at the same time, the solitons inter- 
act with each other. We assume that the pertinent interaction 
potential, Vint, for two (almost black) solitons has the form 
obtained from the exact two-soliton solution of the underlying 
unperturbed NLS system (see details in Ref. l40ll ). namely: 



Vint = 



2sinh [2^/JI(x 2 — x\) 



(21) 



The equation of motion for the two soliton centers (denoted 
hereafter by x\ and x 2 ) can now be obtained from the La- 
grangian: 



1 1 

L = -x\ 2 + -x 2 2 + KffOi) + V e fi(x 2 ) 



(22) 



where the external effective potential V e E is given in Eq. (fT8l) . 
Using the Euler-Lagrange equations, and assuming that the 
two solitons are far away from each other, i.e., \x 2 — x\ | ^> 0, 
we approximate their mutual interaction term in Eq. (l22l) by 
its exponential asymptote and obtain the following equations 
of motion: 



d 2 x\ 
~dfr 



d 2 x 2 
~d& 



-8fi 3 ^e- 2 ^ X2 - Xl) -Gj 2 x x 

e 2 c\ sech 4 (c2£i) tanh(c2#i), (23) 



= S/j J 3/2 e~ 2 ^ X2 ~ Xl) - m 2 



uo x 2 



e 2 c\ sech 4 (c2^2) tanh(c2^2). 



(24) 



where uj = 0/ \/2 is the oscillation frequency of the single 
dark soliton in the absence of W 113 811 . The fixed points X 2 = 
—X\ — X of the above system can be found numerically, 
by setting the left-hand side equal to zero, and employing a 
fixed-point algorithm. 

Next, we can study the stability of these equilibrium posi- 
tions X\ and X 2 , by considering small deviations, r]i and r] 2 , 
around them. By Taylor expanding the nonlinear terms, and 
keeping only linear terms in rji and rj 2 , we derive the follow- 



ing set of linearized equations of motion: 

16/i 2 e- 4 v^ x (77 2 -771) -uj 2 m 



d 2 m 



dt 2 

+ e 2 c x c 2 [5sech 6 (c 2 X) - 4 sech 4 (c 2 X)] rj u (25) 

+ e 2 ac 2 [5sech 6 (c 2 X) - 4 sech 4 (c 2 X)] r/ 2 . (26) 

Let us now consider the normal modes of the linearized sys- 
tem and seek solutions of the form rji = rf^e lUJt where uj is 
the common oscillation frequency of both dark solitons. Then, 
substituting this ansatz into Eqs. d25l)-(l26lh we rewrite the 
equations of motion as a matrix eigenvalue equation, namely: 



,~,2 



■ UJ 2 T] ■ 



/(X) 16/i 2 e- 4 v^ 



16/i 2 e- 4 ^ 



rj. (27) 



where 77 = (7/1, r] 2 ) T and the function f(X) is given by: 

f(X) = -16/i 2 e- 4 ^ x 

+£ 2 cic 2 [5sech 6 (c 2 X) - 4sech 4 (c 2 X)] . (28) 

To this end, it is possible to obtain from Eq. (f27l) the two char- 
acteristic frequencies 



^1,2 



f(X)±W^ 2 e- 4 ^ IX , 



(29) 



where uj\ (uj 2 ) correspond to the in-phase (out-of-phase) os- 
cillations of the two solitons lf38L l4Qh . 

The validity of the above analytical result can now be com- 
pared with the corresponding numerical result obtained by 
the numerical existence theory and the BdG analysis. The 
latter provides the equilibrium configuration, as well as the 
real/imaginary parts of the eigenfrequencies as functions of 
the imaginary potential strength e, as shown in the middle and 
bottom panels of Fig. [3] respectively. In the middle panel, we 
can observe a very good agreement between the positions of 
the theoretically predicted solitary wave centers [cf . the fixed 
points of Eqs. (I23l)-(l24l)l and their numerically exact counter- 
parts. 

For the linearization analysis of the bottom panels, we ob- 
serve the following. In the Hamiltonian case, there exist two 
anomalous modes (and more generally, n anomalous modes 
for a n-soliton state — see, e.g., Ref. [38]). The first, which 
is the lowest nonzero mode in the spectrum, corresponds to 
the in-phase oscillation frequency of the two dark solitons, 
while the second being the fourth finite eigenfrequency for 
our parameters, corresponds to out-of-phase oscillations. The 
two-soliton branch is initially stable up to a critical point (in- 
dicated by the left vertical line), where the first finite mode 
collides with the dipole (alias Kohn) mode [38]; this mode is 
located, for e = 0, at Re(cj) = Q = 0.1. This collision results 
in the emergence of an instability, which is identified by the 
appearance of an imaginary part of the eigenfrequency (see 
the bottom panel of Fig. 0. As e is increased, this mode is 
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eventually detached from the dipole mode and this instabil- 
ity band ceases to exist for a very narrow parametric inter- 
val. However, the quadrupole mode (located, for e = 0, at 
Re(cj) & \/?>VL = 0.17), which has been scattered at about 
£ = 0.08, by the branch just above, now collides with the de- 
tached mode (right vertical line). Upon the latter collision, 
a new (and persistent hereafter) second instability emerges 
(see the bottom panel of Fig. [3]). This way, another quartet of 
complex eigenfrequencies is created, indicating that the two- 
soliton state becomes unstable with a growth rate that contin- 
ues to increase as the gain/loss parameter e increases. 

The analytical result for the two eigenfrequencies ui y2 
[cf. Eq. d29b and solid red and green lines in Fig. [3) essentially 
coincides with the two branches initialized at the two anoma- 
lous modes of the Hamiltonian case, for sufficiently small e; 
however, for larger e, the analytical result is somewhat less ac- 
curate. The observed discrepancy is, at least in part, due to the 
fact that our approximate result is based on a particular ansatz 
of well- separated (thus weakly-interacting), almost black soli- 
tons. However, observing the density profile of the two soliton 
state (top panel of Fig. [3]), one can see that the small density 
dip in the center (which is not included in our ansatz) is cer- 
tainly affecting the interaction between the two solitons; im- 
portantly, this dip becomes larger as e is increased. Thus, for 
relatively large values of e, our perturbative approach may be 
expected to be of lesser value. 

Additionally we note that, as in the case of the single- 
soliton-state, for values e > 0.72 no stationary two-soliton- 
state could be identified. This is again a by-product of the 
nonlinear analogue of the PT-phase transition. More specif- 
ically, the two- soliton branch collides with the three- soliton 
one and they pairwise annihilate, as is occurring to higher (lin- 
ear) eigenvalue pairs in the linear PT-phase transition of [[fl]. 
We illustrate the generality of this effect (and the relevant cas- 
cade of nonlinear eigenvalue collisions) in more detail below. 



B. Three-dark-soliton state 

We now proceed with the investigation of the three- soliton 
branch; examples of the density and phase of such a state are 
shown in the top panel of Fig. HI This state shows similar be- 
havior with the single soliton state, in the sense that the soliton 
located at the center becomes shallower as e is increased. 

The statics and dynamics around equilibria of the three- 
dark-soliton state can be analyzed by the methodology 
adopted in the case of the two-dark- soliton state. Particularly, 
we will determine the equilibrium positions of the three soli- 
tons and study their small- amplitude oscillations around their 
fixed points. First we note that, from symmetry arguments, 
we expect that the fixed points are X 2 = (for the central 
soliton), and X3 = — X\ = X (for each of the two outer soli- 
tons). For this configuration, we can again calculate the fixed 
points for the outer solitons numerically and compare them to 
the particle theory analytical predictions, as illustrated in the 
middle panel of Fig. |U obtaining good agreement for small 
values of e. 

Subsequently, we derive the linearized equations around the 



fixed points, and finally obtain the three relevant eigenfre- 
quencies of the system. This way, we find that the frequencies 
of the normal modes of the system are: 

u>! = y/h(X) (30) 
^2,3 = {\h(X)±f 2 (X)[(h(X)-f 2 (X)f 

i 1/2 -> 1/2 

+8f 3 (X) 2 \ J , (31) 

where functions f%{X) (i = 1,2,3) are given by: 

h(X) = Co 2 -16/i 2 e" 2 v^ 

+ e 2 c lC2 [5sech 6 (c 2 X) - 4 sech 4 (c 2 X)] , (32) 

f 2 (X) = Co 2 - 32/i 2 e- 2 ^ IX + e 2 c x c 2 , (33) 

f 3 (X) = 16fi 2 e- 2 ^ x . (34) 

The above characteristic frequencies will again be compared 
to the eigenfrequencies obtained numerically by means of the 
BdG analysis. For this purpose, in the bottom panels of Fig.H 
we show the BdG linear spectrum of the three-dark-soliton 
branch. 

As observed in the figure, the three-dark-soliton state is ini- 
tially stable (for sufficiently small e), but it becomes unstable 
after a collision of the lowest mode [dashed (red) line] with 
the spectral plane origin, giving rise to the emergence of an 
imaginary pair of eigenvalues. This is directly reminiscent of 
the corresponding instability of the single dark soliton state 
and the mechanism of the instability is expected to persist for 
any configuration with an odd number of dark solitons in this 
system [we have also checked that it arises in the case of a 
five-soliton state]. For larger values of e, another collision be- 
tween the fifth and sixth mode, results in the emergence of 
an imaginary eigenfrequency pair (right vertical dashed line), 
and thus the three soliton state remains unstable. 

At the same time, and as in the single- soliton case, the 
imaginary eigenfrequency corresponding to the first nonzero 
mode, eventually returns to the origin and collides with it, at 
e = 0.72. This is the critical point of the collision of the two- 
soliton and the three- soliton branch and of their pairwise anni- 
hilation. As the relevant critical value of e is approached, the 
center grey soliton of the three- soliton state becomes grayer 
and eventually becomes identical to the central density dip of 
the two- soliton state with which it collides. The above char- 
acteristic frequencies will again be compared to the eigenfre- 
quencies obtained numerically by means of the BdG analysis. 
For this purpose, in the bottom panels of Fig. |4l we show the 
BdG linear spectrum of the three-dark-soliton branch. 

The analytical result for the three eigenfrequencies 
[cf. Eqs. (l30b-d3T1) and solid red, green and magenta lines in 
Fig. SI are in a good agreement with the respective modes 
obtained numerically via the BdG analysis, as long as e is suf- 
ficiently small. In particular, the critical value of e where the 
first finite mode collides with the origin (thus giving rise to the 
onset of the instability) is very well predicted; furthermore, 
the analytically found "trajectory" of the other two branches 
follow, in a fairly good accuracy, the respective numerical re- 
sult, even for moderate and large values of e, with an expected 
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FIG. 4: (Color online) Same as Fig. but for the three soliton state. 
In the bottom panels, solid red, green and magenta lines, correspond 
to the analytical result of Eqs. d30l>-(l371). The value of these modes at 
e = correspond to the first, second and third anomalous modes of 
the Hamiltonian case. 
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FIG. 5: (Color online). The power N as a function of e, for the 
six lowest states. Solid lines correspond to linearly stable states and 
dashed lines to unstable states, with Im(cj) / 0. The vertical lines 
denote the critical values of e found for the single-, two- and three- 
dark- soliton states. This bifurcation diagram illustrates the nonlinear 
PT-phase transition points, at which stationary states disappear in 
pairs through blue- sky bifurcations. 



the bifurcations shown in Fig. [5j is that higher excited states 
sustain these collisions, saddle-center bifurcations and corre- 
sponding blue-sky annihilations for progressively larger val- 
ues of e. This feature motivates us to refer to this process as 
the nonlinear analogue of the "traditional" VT transition, in 
direct correspondence with the pairwise collisions in Ref . Jl|] 
(see, e.g., Fig. 1 of that reference) for the linear setting. Al- 
though the two figures are similar, there are some interesting 
differences. One such is that in the linear picture the lowest 
eigenvalues collide the last, while in the ID nonlinear case the 
lowest eigenstates disappear first 1 50] . This feature has impor- 
tant consequences for the supercritical evolution of the system 
which we will revisit in the dark soliton dynamics section be- 
low. 



progressive discrepancy as e becomes large (when the prob- 
lem is outside the realm of a perturbative treatment). 



C. Nonlinear VT phase transitions 

In order to further elaborate on the emergence of the non- 
linear PT-phase transition, in Fig. [5] we plot the power N as 
function of the parameter e, for the lowest six (ground state 
and the first five excited) states of the system. The top solid 
(blue) branch depicts the stable ground state [i.e., the TF back- 
ground Ub — cf. Eqs. (I9l)-([TT1)1. which ultimately collides with 
the single-dark- soliton branch [dashed (red) line], i.e., the first 
excited state at e « 0.62 (for /i = 3). Furthermore, the next 
pair of branches, corresponding to the two- and three-dark- 
soliton states (depicted by the dashed black and magenta lines, 
respectively), also collide at e « 0.72. This picture of collid- 
ing and annihilating pairs of excited states remains the same 
even for higher excited states — cf. the lowest dashed (green 
and brown) curves in Fig. \5\ showing the annihilation of the 
five- and six-dark- soliton states as well. 

An important general remark, regarding the structure of 



D. The "free space" case 

In our previous considerations, we have studied the case 
where the real part of the external potential was parabolic. 
Nevertheless, and as a complementary setting, it is also inter- 
esting to consider the "free space" case, where the real trap- 
ping potential is absent, i.e., V(x) = 0. Such a situation 
may occur, e.g., in the context of optics where — instead of 
a graded index medium — one may consider a medium with a 
constant linear refractive index, i.e., V(x) = Vq = const.; in 
this case, the pertinent term Vqu that would appear in Eq. (Q]) 
can straightforwardly be removed by a trivial gauge transfor- 
mation. 

An interesting feature characterizing the free space case is 
that stationary multiple dark soliton states can no longer exist: 
in this case, the repulsion between dark solitons cannot be 
counter balanced by any restoring force as the V(x) producing 
the latter is absent. As a result, the most fundamental states 
pertinent to this setting are the ground state and the single- 
dark- soliton state. An approximate solution for the free- space 
ground state uq, can be found employing Eq. (ITTI) and setting 
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FIG. 6: (Color online) Top panel: pertinent eigenfrequency of the lin- 
ear spectrum for the free-space dark soliton with respect to e. Dashed 
(red) line depicts the imaginary eigenfrequency pair arising due to 
the presence of W(x) and the dotted line represents the correspond- 
ing theoretical prediction. Bottom panel: the PT-phase transition 
diagram in the (N, e) plane. Solid (blue) line shows the background 
branch, while dashed (red) line the single-dark- soliton state. The up- 
per and lower insets show characteristic densities and phases for the 
ground state and the dark soliton, respectively. In this case, \i = 1. 



V{x) = 0, namely: 



u 



W 2 exp (2i Wdx 



(35) 



Using the above result, and following the analysis of the pre- 
vious sections, we find that the evolution of the dark soliton 
center is described by [cf. Eq. d20bl: 



cPxo 6 o 

-e xq. 
5 



dt 2 



(36) 



The above equation dictates that a quasi- stationary (almost 
black) dark soliton is always unstable for any finite e, since 
it experiences an effective expulsive force, which tends to dis- 
place it from the origin and set it into motion with an increas- 
ing (magnitude of) acceleration, at least for small values of 
xq. This is in contrast to the case of a parabolic V(x), where 
the respective force exerted on the soliton and the effective ex- 
pulsive force induced by the imaginary potential W{x) could 
balance each other. Thus, the presence of the real parabolic 
trap (or, more generally, of a confining potential) is crucial for 
the existence of stable stationary soliton solutions. 

The top panel of Fig. [6| shows the relevant eigenfrequency 
for the linear stability spectrum of the dark soliton, with re- 
spect to the parameter e. As predicted by Eq. ([36b . as soon 
as e becomes non- vanishing, the dark soliton becomes unsta- 
ble. Its spectrum shows a pair of imaginary eigenfrequen- 
cies emerging from the spectral plane origin in good agree- 
ment (for small/intermediate values of e) with the correspond- 
ing theoretical prediction. Similarly to the case with the real 
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FIG. 7: (Color online) Bifurcation-induced dynamics. Top panels: 
two potential manifestations of the SSB destabilization scenarios for 
a single unstable dark soliton past e = . Bottom panels: soliton 
"sprinklers" spontaneously leading to two or four solitons from the 

(2) 

ground state (initial condition) used for e > Scv . The parameters are 
fj, = 3 and Q = 0.1 and e = 0.4 (top row), e = 0.63 (bottom left), 
and e — 0.66 (bottom right). 



parabolic potential (see Fig. 0, this pair of eigenfrequencies 
eventually returns to — and collides with — the spectral plane 
origin, signaling the termination of the branch and the onset of 
the T^T-transition. In the bottom panel of Fig. [6j the bifurca- 
tion diagram for the ground state and the single-dark- soliton 
in the free-space case is plotted; we observe a behavior similar 
to the one found in the case with the parabolic V(x), but with 
the major difference that beyond the critical point (e = 0.4 
in this case) no stationary states exist, as the only ones such 
have collided and disappeared in the saddle-center bifurcation 
discussed above. 



V. DARK SOLITON DYNAMICS 

We continue by considering the dynamics of the nonlinear 
excited states, by numerically integrating Eq. (Q]) in various 
different regimes in e. To be more specific, our aim is to illus- 
trate the dynamics of the dark soliton upon its destabilization 
after the SSB bifurcation at eil\ and also beyond the nonlin- 
ear PT-phase transition. Generally, to study the dynamics of 
unstable solitons, we prepare a stationary dark soliton solution 
in the unstable regime (dashed parts of the soliton branches in 
Fig. [5]), we add a small random perturbation, and then allow 
the resulting configuration to evolve in time. We have found 
that the dynamics of the system strongly depends on whether 
the soliton will be displaced — as a result of the instability — 
towards x > [the "gain side" of W(#)] or x < [the "loss 
side" of W{x)l 

First, we consider a case where a dark soliton is sponta- 
neously ejected to the gain side of the imaginary potential; a 
pertinent example is shown in the top left panel of Fig. [71 It is 
observed that, in this case, the solitary wave starts to perform 
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oscillations of large amplitude, a behavior which is generic 
for all solitons that are initially located at x > 0. On the other 
hand, if the soliton is initially kicked towards the "lossy" side 
of W(x), upon an initial displacement of a finite ^-dependent 
size, it stops moving and remains quiescent, while the back- 
ground on which it "lives" begins to grow in amplitude and 
width. A relevant example is illustrated in the top right panel 
of Fig. [71 We once again highlight the apparent incompatibil- 
ity of the above scenaria with the expectation of a double well 
effective potential for the soliton and the SSB bifurcation of 
the soliton at xq = 0. This paves the way for the examination 
of the ghost states in the next section. 

The above two cases are the prototypical possibilities (inde- 
pendently of e) for the fate of the dark soliton upon its destabi- 
lization, when the dark soliton state still exists. We now turn to 

(2) 

the examination of what happens beyond e = e C r [where the 
TF background and the single soliton are no longer stationary 
solutions of Eq. (Q}]. There, we use, e.g., an initial condi- 
tion in the form of a TF background u(x,0) = a/// — V(x), 
pertinent to the conservative system 112 811 as a means of explor- 
ing how the system responds when its fundamental states are 
no longer present. We have found that a dark soliton train is 
spontaneously formed, with an increasingly larger number of 
solitons as larger values of e are employed. The bottom pan- 
els of Fig. [7] depict a couple of examples of this phenomenon. 
This constitutes a form of what can be dubbed a "soliton sprin- 
kler". This can, at least in part, also be intuitively connected 
to the observation of Fig. [2 that higher excited multi-soliton 
states persist for larger e than lower ones. Nevertheless, it 
should also be highlighted that the observations typically sug- 
gest that the solitons are nucleated and stay in the vicinity of 
the global minimum of W(x), and more generally tend to be- 
come stationary residing on the "lossy" side of the imaginary 
potential, similarly to the dynamical state observed above and 
accompanied by the same kind of background growth. 

As an additional aspect of the investigation of the dynamics 
of higher-excited (multi-soliton) states, and how their insta- 
bilities manifest themselves when they are subject to small 
perturbations, we consider the runs of the top panel of Fig. [8] 
In the left panel, the dynamics of a two-soliton state is shown. 
At t = the two solitons are in equilibrium, but the soliton on 
the gain side starts to execute oscillations within the trap, as 
was the case before in Fig. [7] for the soliton displaced towards 
the gain side (top left panel of the figure). Thus, inevitably, 
this soliton, upon hitting a turning point in its oscillation, re- 
turns and collides with the other one. The result of the mul- 
tiple collisions that ensue is that one of the solitons continues 
to execute oscillations of increasing amplitude while the other 
one remains effectively stationary on the lossy side. While the 
above sequence of oscillations and collisions is continuously 
repeated, at the same time, the background grows both in am- 
plitude and width. It is interesting to note that the two solitons 
behave in a way that resembles a combination of the two dif- 
ferent scenaria observed in the single soliton case, with the 
relevant dynamics resembling a "superposition" of the ones 
shown in the top two panels of Fig. [7] Yet there is one non- 
trivial difference. While the amplitude of oscillation in the left 
panel of Fig. [7] remains roughly unaltered, the one of the top 
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FIG. 8: (Color online) Bifurcation-induced dynamics. Top panels: 
Manifestation of the dynamical instability of the two soliton state 
(oscillatory instability of the top left panel) and of the three soliton 
state (exponential and oscillatory instability of the top right panel). 
Both end up with one solitary wave sitting on the lossy side and the 
rest impinging on it through collisions between turning points within 
a continuously growing background. The bottom panel shows the 
case without V(x) where the single soliton manifests its instability 
by being set into motion (bottom left panel), while a uniform initial 

(2) 

condition past the Sc r starts emitting dark solitons in the form of a 
sprinkler. The parameters are \i — 3 and Q = 0.1 and e = 0.4 (top 
left), e = 0.4 (top right), while fi = 1 and e = 0.1 (bottom left), and 
e = 0.45 (bottom right). 



left panel of Fig. [8] keeps increasing due to the continuously 
expanding background (which, in turn, shifts the turning point 
of the corresponding oscillation). 

A similar behavior is also found in the case of three soli- 
tons, as shown in the top right panel of Fig. [8j Two of the 
three solitons are now oscillating inside the harmonic poten- 
tial while the third one is trapped in the lossy side of the poten- 
tial. As the solitons pass through the origin, they perform their 
respective oscillations and undergo collisions, in a procedure 
which is continuously repeated; at the same time, once again, 
the background is growing and this, in turn, leads to a growth 
of the turning points of the nearly alternating oscillations of 
the two non- stationary solitons. 

We now compare the above results to the ones correspond- 
ing to the free-space case, i.e., when the trapping potential 
V(x) is absent. Recall that in this case only two stationary so- 
lutions exist (the ground state and the single soliton) and the 
soliton is always unstable. Since the background wave is now 
of constant density, when the soliton is perturbed, it will start 
to move away from the origin and will accelerate according to 
Eq. ([36b up to the point where it is practically no longer under 
the influence of the (exponentially localized) imaginary part 
of the potential. Thereafter, it will acquire a nearly constant 
velocity and will keep moving in that direction, as shown in 
the bottom left panel of Fig. [8] (the absence of the real part of 
the potential leads to a lack of a turning point in this case). In 
the absence of V(x), the soliton dynamics does not depend on 
which direction (towards x < or x > 0) the soliton will be 
initially displaced. An initial kick in the direction of x < 
will produce a similar phenomenology but for negative values 
of x. 

Finally, we consider the parameter regime e > e cr , i.e., 
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after the occurrence of the PT-phase transition where no sta- 
tionary solutions exist. We numerically integrate Eq. (Q]), with 
an initial condition of the form of a plane wave \u(x, 0)| = 
y/JI, and the result is shown in the right bottom panel of Fig. [8] 
The initial plane wave begins to emit dark solitons, which are 
all starting to propagate along the "lossy" direction. Contrary 
to the confined case where the trapping potential was keeping 
the solitons near the origin, the solitons are continuously cre- 
ated and are now free to propagate along the x axis. Hence, 
once the "sprinkling" process begins, it continues indefinitely. 

VI. GHOST STATES AND THEIR DYNAMICAL ROLE 

It is now relevant to try to address some of the remaining 
questions, concerning the nature of the daughter states within 
the pitchfork bifurcation (of Fig. [2] and by extension e.g. in 
Fig. |4]) and the unexpected asymmetric behavior of the dy- 
namical evolution of unstable dark solitary waves in Fig. [71 
A related issue concerns the feature that solitary waves be- 
come stationary on the lossy side of the potential (while the 
corresponding background starts growing); see Fig. [7] and by 
extension also Fig.[8j 

A first way to try to address this problem is to consider 
the possibility (as suggested by the direct numerical simula- 
tions) of a state involving a stationary dark solitary wave off 
of xo = (and, in particular, for xq < 0). Then the steady 
state problem would read: 

\±u = Cu + iW(x)u + \u\ 2 u, (37) 

where C is used to denote the linear, Hermitian part of the 
right hand side operator. The corresponding conjugate equa- 
tion then reads: 

H*u* = Cu* - iW(x)u* + M V. (38) 

Now, by multiplying Eq. (l37l) by u*, and Eq. (l38l) by u, inte- 
grating and subtracting the second equation from the first, we 
obtain the self-consistency condition for the imaginary part of 
the chemical potential: 

Im(/i) = J W(x)\u\ 2 dx. (39) 

Interestingly, the right hand side is directly related to the one 
arising in Eq. © and, thus, also controls the time dependence 
of the atom number/optical power (or mathematically the L 2 
norm). 

The key realization that emerges from Eq. d39b is that since 
\u\ 2 is an even function for solitary structures centered at 
xo = (and W(x) is anti- symmetric), the only way for our 
states to be centered at xo ^ is if /ij = lm(fi) ^ 0. How- 
ever, this is an atypical feature for states arising in the context 
of the nonlinear Schrodinger (NLS) equation. More specifi- 
cally, a fundamental underlying assumption behind the iden- 
tification of stationary states, critically employing the U(l) 
invariance of the model, is the use of the product ansatz for 
stationary states of the form u(x,t) = e~ ZfJjt f(x). Now, we 
find ourselves at a junction where solutions may exist with 



a complex fi. This, more specifically, implies that they will 
be indeed exact solutions of the stationary NLS problem, yet 
their incompatibility with the above ansatz precludes them 
from being exact solutions of the original dynamical problem 
ofEq. CD). 

The natural next questions then involve whether these states 
can, in fact, be identified as such (exact solutions of the steady 
state problem) and, perhaps more importantly, what is their 
role in the observed dynamics. The first of these questions is 
answered in Fig. [9J where we have performed a continuation 
of these states as a function of e, starting from the Hamilto- 
nian limit in the absence of gain/loss. The top left panel of 
the figure presents the center of the solution as a function of 
the parameter, while the top right one evaluates the imaginary 
part of the chemical potential, by self-consistently enforcing 
the condition ([39b . It is key to mention here that for the so- 
lutions with negative xo, it is straightforward to see that they 
correspond to a setting with /ij > [from Eq. d39bl. On 
the other hand, the term e~ Zflt will then involve a part associ- 
ated with e - *( z/iJ ) £ and hence this state will be associated with 
growth over time. On the other hand, the state with /ij < 
will be connected to decay over time and does not appear to 
be relevant for physical observations. It is for that reason that 
one of the two states (the robust one) is denoted with a bold 
solid line in the top panels of Fig. [9l while the other one is 
denoted by a thin line (and unstable states centered at xq = 
in the diagram are denoted by a dashed line). 

First of all, this observation restores the canonical (and ex- 
pected) pitchfork nature of the observed bifurcation (see top 
panels of Fig. [9] and the daughter state of the bottom left 
panel). But, at the same time it reconciles that nature with 
the absence of any "regular" daughter states. A stability anal- 
ysis of these states is not particularly meaningful per se (as 
these are not stationary states of the original problem), yet it 
is natural to expect that it will be suggestive regarding the evo- 
lution of perturbations in the vicinity of such states. This, as 
shown in the bottom right panel for the displaced state of the 
bottom left one in Fig. illustrates the effective stability of 
the resulting ghost state. 

We now turn to the dynamical implications of the existence 
of such ghost modes. The asymmetry between the growth of 
the mode with the dark soliton at xo < and the decay of 
the one with xq > justifies the observed asymmetry of the 
dynamical evolution of the perturbations in Fig. [7] To illus- 
trate the relevance of the ghost states in the case where the 
random perturbation manifests the instability of the xo = 
soliton by kicking it to the left, we performed some relevant 
numerical experiments in Fig.[l0l for e = 0.15. Specifically, 
by initiating the soliton at xq = 0, we confirm that it goes and 
sits precisely at the location (for that particular e) where the 
ghost state has the soliton centered (cf. with the out of scale 
profile of the shifted soliton within the graph). On the other 
hand, for a clear comparison, in the right panel of Fig.[l0l we 
have evolved the dynamical equations with exactly that proper 
ghost state for the particular e as initial condition (wherein the 
soliton is located at —1.5 for e = 0.15; see the right panel of 
the figure). The dynamics clearly illustrates that the soliton 
stays immobile (as expected by the nature of the ghost state), 
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FIG. 9: (Color online) Properties of the ghost states. The top left 
panel illustrates their (off-center) location with the bold line illus- 
trating the dynamically robust one. The center line corresponds to 
the dark soliton at xo = and the solid to dashed transition cor- 
responds to the SSB bifurcation. Another diagram illustrating the 
pitchfork character of the transition is given in the top right panel 
through the imaginary part of the chemical potential /xj as a func- 
tion of the gain/loss parameter. Again, the persisting state is given 
in bold. A typical example of the squared density and phase (solid 
and dashed, respectively) of a ghost state for e = 0.15 is given in 
the bottom left panel. Finally, a suggestive linear stability analysis 
(as if this were a true steady state of the system) is performed on the 
bottom right suggesting also the dynamical robustness of the relevant 
state due to the absence of eigenvalues with a positive real part (i.e., 
of ones corresponding to growth). 



while the background increases in amplitude and also the 
width of the solution increases. These latter growth features 
are strongly reminiscent of the evolution of the left panel of 
the figure past the point of destabilization and convergence to 
a solitary wave centered at xo = — 1.5. To further cement this 
proximity, we have computed the growth of the background 
amplitude i.e., the maximum absolute value of the field A(t) 
and have directly compared it to the simplistic law (implied by 
the ghost state growth rate ///) A(t) = A(0) exp(/ijt), yield- 
ing the kind of agreement shown in the bottom panel of the 
figure. The natural conclusion is that while the ghost states 
are not genuine (steady) attractors of the original dynamics, 
nevertheless the evolution of the system closely shadows the 
relevant states and utilizes their asymmetry. 

As for the case where the dark soliton moves to xq > 
due to the perturbation, then the ghost state is simply unten- 
able as a dynamical state (due to its decay in time), hence the 
dark soliton proceeds to escape from the region of influence 
of the gain and moves towards a turning point. This motion 
imparts sufficient momentum to the solitary wave that it can 
subsequently execute oscillations around xo = 0, as shown in 
Fig. [7] (top right). Moreover, dynamical features such as the 
generation of multi-solitons sitting on the lossy side can be ex- 
plained upon the consideration of ghost states associated with 
multi-soliton states (see e.g. the 3-soliton state of the bottom 
right of Fig. 0, while features such as those of the top panel 
of Fig. [8] can be explained as well. The latter can be compre- 
hended on the basis of a superposition of a ghost state centered 
on the lossy side with a growing background, and of one (as 
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FIG. 10: (Color online) The top left panel shows the space-time evo- 
lution of the contour plot of the unstable dark soliton at xo — 
for e — 0.15. For comparison (and out of scale), we plot (using a 
solid [blue] line) the exact obtained ghost state for the same gain/loss 
parameter. Clearly the dynamics gets attracted to this state, as is 
also corroborated by its evolution as initial condition in the top right 
panel (the dark soliton stays put, but the background width and am- 
plitude grows in a way directly analogous to the case on the left). 
Finally, the comparison of the amplitude background growth A(t) 
(solid line) with the law implied by the growth rate of the ghost state 
A(t) = A(0) exp(/i/t) (dashed) is shown on the bottom panel, sug- 
gesting good agreement between the two. 

in top left, or more as in the top right panel of Fig. [TO)) soli- 
ton(s) moving to the opposite direction and hence executing 
oscillations due to the absence of a ghost state to trap it. The 
above characteristics of the ghost states thus naturally lend 
themselves to the explanation of the static and dynamic phe- 
nomenology observed herein. 

VII. SYMMETRY BREAKING AND NONLINEAR VT 
PHASE TRANSITIONS IN TWO DIMENSIONS. 

In order to consider the generalization of the ideas pre- 
sented herein to higher dimensions, we now consider the case 
of a 2D VT- symmetric potential [W(—x, —y) = —W(x,y)] 
with 

W(x, y) = e(x + y)e~ {x2+y2)/A . (40) 

The configurations sustained in the 2D setting follow a similar 
symmetry breaking bifurcation scenario as in their ID coun- 
terparts. However, the bifurcation structure is naturally ex- 
pected to be more complex for the higher dimensional case. 
Figure [TT] displays the bifurcation scenario for 2D solutions 
with and without topological charge for Vt = 0.2. The bi- 
furcation diagram includes the chargeless ground state (the 
TF background cloud), solutions bearing from one to six vor- 
tices, and the dark soliton stripe. Similar to the ID case, the 
TF background is stable in all its domain of existence and col- 
lides, as £ increases, with an excited state (the two- vortex so- 

(2) 

lution or vortex-dipole) in a blue-sky bifurcation at e = £cr . 
However, in contrast to the ID case where this collision oc- 
curs with the first excited state, in the 2D case the collision 
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FIG. 1 1 : (Color online) Bifurcation diagram of stationary profiles in 
the 2D case. Stable and unstable branches, as determined by the cor- 
responding BdG analysis, are depicted, respectively, with solid and 
dashed lines. The large symbols (circles, squares and triangles) de- 
note the location in parameter space for the series of profiles depicted 
in Fig. Q2] The parameter values are: — 2 and Q = 0.2. 



happens with the second excited state consisting of a vortex- 
dipole (a vortex pair with opposite charges). This is due to the 
fact that the TF background has no topological charge and, in 
turn, the vortex-dipole has also no net topological charge. This 
allows for the emergence of the vortex-dipole from a central 
dip on the TF background [see top panel of Fig. 02 a)] that 
becomes deeper as e increases. This blue- sky bifurcation is 
depicted in Fig. [TT] where the stable (center) TF branch (solid 
blue curve denoted by TF) collides with the unstable (saddle) 
vortex-dipole branch (red dashed curve denoted by2-vort) 

(2) 

at the ID equivalent of the critical point e = £cr > 

It is possible to further follow the vortex-dipole branch by 
decreasing e where we note that it regains stability for suffi- 
ciently smaller values of e (see transition between the dashed 
and solid portions of the red line for the branch denoted by 
2-vort). As this branch is followed further (from top to 
bottom in Fig.lTTb a series of bifurcations arise where the vor- 
tices are drawn towards the edge of the cloud, a central dip 
deepens leading eventually to the emergence of a new vor- 
tex dipole in the middle of the cloud (namely, the emergence 
of a higher excited state). Therefore, the branches with even 
number of vortices are all inter-connected in this bifurcation 
scenario and only the TF branch and part of the vortex-dipole 
branch are stable. However, as more and more vortex pairs 
emerge, the extent of the cloud "saturates" as it can no longer 
support new vortex pairs and therefore is replaced by a dark 
soliton stripe (see dashed blue line denoted by DS in Fig.fTTT). 
This overall bifurcating scenario for even number of vortices 
is symmetric such that if £ — ^ — £ the solutions are just flipped 
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FIG. 12: (Color online) Stationary states in the 2D case correspond- 
ing to the main bifurcation branches in Fig. [TT] The respective left 
(right) columns depict the density (phase) profiles of the solutions, 
(a) Profiles for the branch containing only even number of vortices 
corresponding to the circles in Fig. [TT] [from top to bottom], (b) Pro- 
files for the branch starting with a single vortex and connecting with 
three vortex branch corresponding to the squares in Fig. [TT] [from top 
to bottom], (c) Profiles for the branch starting with three symmetric 
vortices and ending with four vortices corresponding to the triangles 
in Fig. ITTI [from top to bottom]. The parameter values are: \i = 2 
and Q = 0.2. 



by (x,y) —> {—x,—y). Figure [T2la) depicts the density and 
phase profiles for different representatives along this branch. 

Let us now describe the bifurcation scenario for odd num- 
ber of vortices. The first excited state, corresponding to a sin- 
gle vortex, starts with a single vortex at the origin for e = 
and it is stable for small values of e. As in the ID case, this 
first excited state suffers a SSB bifurcation at £ — £i) where 
it loses its stability. As e is increased, the single vortex moves 
towards the periphery of the cloud and a dip at the center of 
the cloud deepens until a vortex dipole emerges at the center in 
the same way as new vortex-dipoles emerged for the branches 
with even number of vortices described above. This scenario 
connects, again through a blue-sky bifurcation, the one- vortex 
branch (magenta line denoted by 1-vort in Fig. ITTb with 

the asymmetric three- vortex (H h vortex tripole) branch 

(green dashed line denoted by 3-vort in Fig.Qj]). The se- 
ries of snapshots at the parameters depicted by the squares 
in Fig. [TT]is depicted in Fig. [T2lb). It is nevertheless rele- 
vant to note that the blue- sky bifurcation for the single-triple 
vortex state (and more generally for any higher order pair of 
vortices) happens for values of £ lower than that of the bifurca- 
tion involving the TF ground state and the vortex-dipole state. 
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FIG. 13: (Color online) Dynamics of unstable steady states in the 2D 
case. All panels depict the density profiles at the times indicated for 
\i — 2 and Q = 0.2. The thin diagonal line in all panels depicts the 
border between the loss and gain sides of W(x,y). (a-b) Unstable 
single vortex state for £ = 0.2 and e = 0.35, respectively, (c) Stable 
single vortex configuration for e = 0.1 perturbed by a small rota- 
tion of 5 degrees, (d) Stable single vortex configuration for £ = 0.1 
perturbed by a small rotation of 4 degrees performs a periodic oscil- 
lation around its original position (a single period is shown, see zoom 
in the inset). 



Interestingly, although the structure of the relevant bifurca- 
tions is reminiscent of its one-dimensional analogue, there are 
also non-trivial differences, including the topological charge 
issue mentioned above, as well as the order in which the bi- 
furcations occur (between more fundamental and more highly 
excited states) which is reversed. As it is evident from the fig- 
ure, the asymmetric three- vortex branch eventually connects 
with the symmetric one for values of e — )> 0. A similar bifur- 
cation occurs with the symmetric three- vortex branch, which 
becomes asymmetric with a deepening dip at the center where 
a vortex pair emerges (at the same time that a vortex is lost 
at the periphery), connecting in this way with the four- vortex 
branch [see series of snapshots in panels H2c)]. As a relevant 
aside, let us also mention that all of the above true stationary 
states of the system (with real chemical potential or propaga- 
tion constant) clearly form along the diagonal x + y = on 
which there is no gain or loss according to the prescription of 



FIG. 14: (Color online) Same as in Fig. [T3l for the following two- 
vortex scenarios, (a) Unstable tight vortex-dipole state belonging to 
the branch directly connecting with the TF cloud (depicted by the 
red dashed line denoted by 2-vort in Fig. [TT) for e = 0.4. (b) 
Unstable, well separated, two-vortex state belonging to the branch 
connecting with the four- vortex branch (depicted by the black dashed 
line denoted by 2-vort in Fig.Qj} for e — 0.2. 



It is important to mention that the precise structure of the bi- 
furcation diagram depends on the values of the trap strength Vt 
and the propagation constant fi. In general terms, for weaker 
Q and/or larger fi the spatial extent of the TF background will 
be larger and thus allowing for a longer bifurcating chain of 
higher-order vortex states (before all the vortices finally merge 
into a dark soliton stripe). Nonetheless, the displayed SSB in- 
stabilities and the nonlinear VT phase transition involving the 
cascade of blue-sky bifurcations appear to be universal in con- 
fining T^T-symmetric potentials as considered herein. 

Let us now describe the dynamics of the unstable steady 
states described in the above bifurcation scenario. We first 
start by describing the unstable dynamics of single vortex 
states. In this case, in analogy with the ID case, single vortex 
states tend to migrate towards the lossy side of the potential. 
This tendency is depicted by the series of density snapshots 
in Figs. H2a) and fT3lb) where unstable single vortex states 
for e — 0.2 and e = 0.35, respectively, migrate to differ- 
ent positions within the lossy side. This phenomenology can 
be very clearly understood on the basis of our ghost state in- 
terpretation. Once again, ghost states are present here corre- 
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FIG. 15: (Color online) Same as in Fig. [13] for the following three- 
vortex scenarios, (a-b) Unstable asymmetric three- vortex states (de- 
picted by the green dashed line denoted by 3-vort in Fig. [TTJ for 
e — 0.2 and e — 0.35, respectively, (c) Unstable symmetric three- 
vortex state for £ = 0. (d) Unstable asymmetric three- vortex state 
(depicted by the green dashed line denoted by 3-vort in Fig. [TTJ 
fore = 0.1. 



sponding to a unique shift of the vortex center to a location 
such that xo + yo < 0. This state is associated with /ij > 0, 
due to the fact that J W(x, y)\u(x, y)\ 2 dxdy > and hence 
leads to growth of the amplitude and width of the condensate, 
a trait which can be clearly observed in the relevant panels. 
It is also important to mention that the stable single vortex, 
close to where it loses stability, is only weakly stable and can 
be 'kicked out' by a relatively small perturbation towards the 
lossy side. This scenario is depicted in Fig. Etc) where the 
stable single vortex configuration for e = 0.1 is perturbed by 
applying a small rotation of 5 degrees. Although the rotation 
is very small and that the unperturbed configuration is stable, 
we note that the vortex eventually migrates (after a long tran- 
sient, see that time in the series of snapshots runs to t=3000) 
towards the lossy side. Further numerics shows that angles 
larger than 4 degrees display the same dynamics while smaller 
angles correspond to the single vortex returning periodically 
towards its original position [see for example the periodic or- 
bit depicted in Fig. HHd) where the stable vortex is initially 
rotated by 4 degrees] . 

Let us now describe the evolution of unstable vortex- 



FIG. 16: (Color online) Same as in Fig.[l3]for the following scenar- 
ios, (a-b) Unstable dark soliton stripes (depicted by the lowest blue 
dashed line denoted by DS in Fig. [TTJ for e = 0.1 and e = 0, respec- 
tively, (c) Thomas-Fermi (ground) state obtained for s = 0.4 and let 
to evolve for e — 0.43. 



dipoles. In Fig.[l4ta) we depict the evolution of a tight vortex- 
dipole close to the parameter values where it is created (i.e., 
in the dip at the center of the cloud). As the figure shows, 
the tight vortex-dipole undergoes a small excursion into the 
lossy side and, as time progresses, it re-approaches to the cen- 
ter of the cloud at the same time that the cloud grows (as in 
the ID case). Clearly the relevant interpretation here involves 
the existence of a ghost state which attracts the dynamics and 
contains both vortex members of the dipole on the lossy side 
of the system. In Fig. [Bib) we depict the dynamics arising 
from a well separated unstable vortex-dipole state. We typi- 
cally observe that, for a well separated unstable vortex dipole, 
one vortex is trapped at a suitable location within the lossy 
side while the other vortex circles around slowly approaching 
the periphery of the cloud where it is finally absorbed. This 
is in a number of ways tantamount to the two- soliton dynam- 
ics that we explored in Fig. [8] Generally, it is clear that the 
3 possible dynamical scenaria for the vortices constituting the 
dipole and involving ghost states here are that either both vor- 
tices transition to the loss side rFig.PPfta)!, or one rFig.[T4Tb)l 
or none (not shown here). 

Let us now briefly describe some typical evolution exam- 
ples for higher order unstable states. In particular, let us start 
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by discussing the three- vortex states. Figure [T5[a) depicts the 
evolution of an unstable asymmetric 3 -vortex profile. In this 
case, the outer vortex migrates towards the attracting basin 
within the lossy side while the remaining vortex-dipole gets 
violently ejected and the two vortices annhilate each other. 
This again confirms the dynamical relevance of the single- 
vortex corresponding ghost state. Another possible scenario 
for the 3 -vortex configuration is depicted in Fig.[T3b) where 
this time the central vortex-dipole slowly migrates towards the 
lossy side and slowly re-approaches the center of the trap as 
time progresses while the outer vortex performs large excur- 
sions around the periphery of the cloud where it is eventu- 
ally absorbed. In this setting, we observe once again the dy- 
namics of a vortex-dipole ghost state (with the background in 
this, as well as in the previous example featuring the analyzed 
growth and concurrent spreading). For contrasting purposes, 
we depict in Fig.[l5tc) the dynamics of an unstable symmetric 
three- vortex state for e — 0. In this case, the external po- 
tential is conservative and the motion of the three vortices is 
Hamiltonian and thus, typically, the three vortices will remain 
orbiting inside the cloud for very long times. Nevertheless, 
notice in the latter case the well established (and generic) in- 
stability of the Hamiltonian vortex tripole which leads to its 
break-up into a vortex dipole and a single vortex 00. Fi- 
nally, in Fig.[l5td) we depict the evolution of another unstable 
asymmetric state. In this case it can be observed that two vor- 
tices annhilate each other and the third does not get absorbed 
by the lossy side as in Fig. [I5ta), but it performs large oscil- 
lations close to the periphery of the cloud where it is finally 
absorbed. The precession of this vortex due to its initial place- 
ment within the gain side suggests a direct analogy of this case 
(and after the dipole annihilation) with the top left panel of 
Fig.13 

Finally, we briefly explore the dynamics of unstable dark 
soliton stripes. Dark soliton stripes are known to be unstable 
in the Hamiltonian case (e = 0) due to the so-called snaking 
instability H^, |47|] (although they can be rendered stable by 
a very tightly confining trap 114 811 or a wall-like external po- 
tential |49|]). Figure H3a) depicts the typical evolution of an 
unstable dark soliton stripe solution in the non-Hamiltonian 
(PT-symmetric) case. In this case, the dark soliton stripe de- 
cays into a chain of alternating charge vortices induced by 
the snaking instability and their long term behavior includes 
annihilation and large oscillations where they get eventually 
absorbed at the periphery of the cloud leaving behind a highly 
perturbed (chargeless) Thomas-Fermi ground state. In con- 
trast, as depicted in Fig. [T6lb), for the Hamiltonian case of 
e = 0, most of the vortices that formed from the snaking in- 
stability remain interacting for long times (i.e., the latter are 
more likely to survive here, analogously to what we also saw 
for the tripole case above). Finally, in Fig.QjJc) we depict the 
evolution when starting with a TF cloud for e past the critical 
value of existence of the stationary TF cloud. In this example, 
we compute the stationary (stable), chargeless, steady TF pro- 
file for e = 0.4 and run it with e = 0.43 (i.e., past the blue- sky 
bifurcation between the TF cloud and the vortex dipole states). 
We observe in this simulation that the central dip of the TF 
(which contains no charge) rapidly gets converted into a tight 



vortex-dipole state which does a small excursion towards the 
lossy side [very similar to the dynamics of the unstable tight 
vortex-dipole profile of Fig.[l4ta)], clearly suggesting the per- 
sistence and dynamical selection of the corresponding vortex- 
dipole ghost state. This behavior is akin to what we observed 

(2) 

in ID for values of e past £cr (the dark soliton 'sprinkler' 
case) where excited states with topological charge are dynam- 
ically created from chargeless configurations. Hence, the two- 
dimensional setting may be used as a "vortex sprinkler". 



VIII. CONCLUSIONS. 

In the present work, we offered a basic, yet self-contained 
view on the fundamental modifications that the phenomenol- 
ogy of nonlinear entities sustains in the presence of VT- 
symmetric confining potentials (although the case where there 
is no real potential was also touched upon as well). The pro- 
totypical among the observed features of such systems were 
found to be the emergence of symmetry breaking bifurcations 
that destabilizes single, triple, five, etc. soliton states in one 
dimension and similarly one- or two- or higher vortex states in 
two-dimensions. The second important and unexpected fea- 
ture was the emergence of a nonlinear analogue of the VT- 
phase-transition, whereby e.g. the ground and first excited, the 
second and third excited, the fourth and fifth excited state and 
so on (in ID) are led to pairwise collisions and blue- sky bi- 
furcations. Generalizations of these effects which, however, 
properly respect the topological nature of the states (e.g. the 
zero-vortex with the vortex dipole, the single vortex with the 
vortex tripole etc.) are found to arise in the two-dimensional 
setting. In addition, the pitchfork nature of the symmetry- 
breaking bifurcation was illustrated through the examination 
of the so-called ghost states and the proof, not only of their 
genuinely complex eigenvalue parameter, but also of their di- 
rect relevance for the dynamics of the full problem, despite 
the fact that they do not constitute exact solutions thereof (but 
rather only solutions of the static problem). The dynamical 
evolution of unstable soliton and vortex states was examined 
for all the identified instabilities, revealing a wealth of pos- 
sibilities involving ghost states, soliton/vortex sprinklers, and 
oscillations between the loss and gain sides. 

It would be interesting to offer a detailed map of the ghost 
states of the system. Here, we have illustrated their existence 
and observed many of their dynamical implications, yet a sys- 
tematic characterization of their emergence (and perhaps the 
more challenging understanding of their basins of attraction) 
certainly merits further examination. On the other hand, a 
feature that has been examined in some detail in the context 
of PT-symmetric linear problems which is, arguably, worth- 
while to consider here concerns the analytic continuation of 
the states past the point of the PT-phase-transition. This was 
already touched upon in the nonlinear realm for double well 
potentials in the work of Refs. Il32l - l35n . but it would be rele- 
vant to offer a definitive view of such analytically continued 
states for general potentials. Finally, a canonical set of in- 
vestigations which is still missing concerns the effects of such 
potentials in three-dimensional continuum or one- and higher- 
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dimensional (infinite) lattice contexts. These topics will be pursued in future publications. 
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